A Classical Density-Functional Theory for Describing Water Interfaces 



Jessica Hughes, Eric Krebs, and David Roundy 

Department of Physics, Oregon State University, Corvallis, OR 97331, 



USA 



(N 

o 
3 



o 

c3 



o 
o 



> 
in 
m 
O 
q 

ON' 

o: 

(N 



X 



We develop a classical density functional for water which combines the White Bear fundamental- 
measure theory (FMT) functional for the hard sphere fluid with attractive interactions based on 
the Statistical Associating Fluid Theory (SAFT-VR). This functional reproduces the properties of 
water at both long and short length scales over a wide range of temperatures, and is computationally 
efficient, comparable to the cost of FMT itself. We demonstrate our functional by applying it to 
systems composed of two hard rods, four hard rods arranged in a square and hard spheres in water. 



I. INTRODUCTION 

A large fraction of interesting chemistry — including all 
of molecular biology — takes place in aqueous solution. 
However, while quantum chemistry enables us to calcu- 
late the ground state energies of large molecules in vac- 
uum, prediction of the free energy of even the smallest 
molecules in the presence of a solvent poses a continuing 
challenge due to the complex structure of a liquid and 
the computational cost of ab initio molecular dynam- 
ics^. The current state-of-the art in ab initio molecular 
dynamics is limited to a few hundred water molecules per 
unit celli On top of this, standard ab initio methods us- 
ing classical molecular dynamics without van der Waals 
corrections strongly over-structure water, to the point 
that ice melts at over 400K— ! There has been a flurry 
of recent publications implicating van der Waals effects 
as significant in reducing this over-structuring^—. It has 
also been found that the inclusion of nuclear quantum 
effects can provide similar improvements^. Each of these 
corrections imposes an additional computational burden 
on an approach that is already feasible for only a very 
small number of water molecules. A more efficient ap- 
proach is needed in order to study nanoscale and larger 
solutes. 



A. Classical density-functional theory 

Numerous approaches have been developed to approx- 
imate the effect of water as a solvent in atomistic cal- 
culations. Each of these approaches gives an adequate 
description of some aspect of interactions with water, 
but none of them is adequate for describing all these 
interactions with an accuracy close to that attained by 
ab initio calculations. The theory of Lum, Chandler and 
Weeks (LCW)i£, for instance, can accurately describe the 
free energy cost of creating a cavity by placing a solute 
in water, but does not lend itself to extensions treating 
the strong interaction of water with hydrophilic solutes. 
Treatment of water as a continuum dielectric with a cav- 
ity surrounding each solute can give accurate predictions 
for the energy of solvation of ionsii~— , but provides no 
information about the size of this cavity. In a physically 
consistent approach, the size of the cavity will naturally 
arise from a balance between the free energy required to 



create the cavity, the attraction between the water and 
the solute, and the steric repulsion which opens up the 
cavity in the first place. 

One promising approach for an efficient continuum de- 
scription of water is that of classical density-functional 
theory (DFT), which is an approach for evaluating the 
free energy and thermally averaged density of fluids in 
an arbitrary external potentiate. The foundation of 
classical DFT is the Mermin theorem^, which extends 
the Hohcnberg-Kohn theorem^ to non-zero temperature, 
stating that 



A(T) = min 
n(r) 



F[n(r),T}+ / V ext {v)n(v)dr 



(1) 



where A(T) is the Hclmholtz free energy of a system in 
the external potential V ex t at temperature T, n(r) is the 
density of atoms or molecules, and F[n(r),T] is a uni- 
versal free-energy functional for the fluid, which is inde- 
pendent of the external potential V ex p Classical DFT 
is a natural framework for creating a more flexible the- 
ory of hydrophobicity that can readily describe interac- 
tion of water with arbitrary external potentials — such as 
potentials describing strong interactions with solutes or 
surfaces. 

A number of exact properties arc easily achieved in 
the density-functional framework, such as the contact- 
value theorem, which ensures a correct excess chemical 
potential for small hard solutes. Much of the research 
on classical density-functional theory has focused on the 
hard-sphere fluid22r— , which has led to a number of so- 
phisticated functionals, such as the fundamental- measure 
theory (FMT) functionals^—. These functionals are en- 
tirely expressed as an integral of local functions of a few 
convolutions of the density (fundamental measures) that 
can be efficiently computed. We will use the White Bear 
version of the FMT functional^. This functional reduces 
to the Carnahan-Starling equation of state in the homo- 
geneous limit, and it reproduces the exact free energy in 
the strongly-confined limit of a small cavity. 

A number of classical density functionals have been 
developed for water— ~— , each of which captures some of 
the qualitative behavior of water. However, each of these 
functionals also fail to capture some of water's unique 
properties. For instance, the functional of Lischner et 
aJ2& treats the surface tension correctly, but can only be 
used at room temperature, and thus captures none of 
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the temperature-dependence of water. A functional by 
Chuev and SkoloA* 3 - 5 - uses an ad hoc modification of FMT 
that can predict hydrophobic hydration near tempera- 
tures of 298 K, but does not produce a correct equation 
of state due to their method producting a high value 
for pressure. A number of classical density function- 
als have recently been produced that are based on Sta- 
tistical Associating Fluid Theory (SAFT^r^MLSE-^. 
These functionals are based on a perturbative thermo- 
dynamic expansion, and do reproduce the temperature- 
dependence of water's properties. 

B. Statistical associating fluid theory 

Statistical Associating Fluid Theory (SAFT) is a the- 
ory describing complex fluids in which hydrogen bonding 
plays a significant role 41 -. SAFT is used to accurately 
model the equations of state of both pure fluids and 
mixtures over a wide range of temperatures and pres- 
sures. SAFT is based on Wertheim's first-order thermo- 
dynamic perturbation theory (TPTl)^s— , which allows 
it to account for strong associative interactions between 
molecules. 

The SAFT Hclmholtz free energy is composed of five 
terms: 

F Fid ~\~ Ffis ~f" Fdisp T F assoc -\- Fchain^ (2) 

where the first three terms — ideal gas, hard-sphere repul- 
sion and dispersion — encompass the monomer contribu- 
tion to the free energy, the fourth is the association free 
energy, describing hydrogen bonds, and the final term is 
the chain formation energy for fluids that are chain poly- 
mers. While a number of formulations of SAFT have 
been published, we will focus on SAFT-VB3, which was 
used by Clark et al to construct an optimal SAFT model 
for water—. All but one of the six empirical parameters 
used in the functional introduced in this paper are taken 
directly from this Clark et al paper. As an example of the 
power of this model, it predicts an enthalpy of vaporiza- 
tion at 100°C of AH vap = 39.41 kJ/mol, compared with 
the experimental value AH vap = 40.65 kJ/mol 4 ^, with an 
error of only a few percent. We show a phase diagram for 
this optimal SAFT model for water in Figure [TJ which 
demonstrates that its vapor pressure as a function of tem- 
perature is very accurate, while the liquid density shows 
larger discrepancies. The critical point is very poorly de- 
scribed, which is a common failing of models that are 
based on a mean-field expansion. 

SAFT has been used to construct classical density 
functionals, which are often used to study the surface 
tension as a function of temperature 2 ^ 34 ' 37-39 ' 42 ' 43 . Such 
functionals have qualitatively predicted the dependence 
of surface tension on temperature, but they also overesti- 
mate the surface tension by about 50%, and most SAFT- 
bascd functionals are unsuitcd for studying systems that 
have density variations on a molecular length scale due to 
the use of a local density approximation 3 ^ 32 ' 34 ' 38 ' 39 ' 42 ' 43 . 
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FIG. 1. The pressure versus density for various temperatures, 
including experimental pressure data from NIST—. The solid 
colored lines indicate the computationally calculated pressure 
and the dotted colored lines are NIST data points. The solid 
and dotted black lines represent the theoretical and experi- 
mental coexistence curves. 



Functionals constructed using a local density approx- 
imation fail to satisfy the contact-value theorem, and 
therefore incorrectly model small hard solutes. The 
contact-value theorem states that the pressure any fluid 
exerts on a hard wall interface is proportional to the con- 
tact density of the fluid 50 : 

p = n c k B T, (3) 

where n c is the contact density, ks is the Boltzmann 
constant and T is the temperature of the fluid. This leads 
to the property that the excess chemical potential of a 
small hard solute is proportional to the solvent-excluded 
volume: 

F = nk B TV. (4) 

The contact-value theorem is satisfied by classical density 
functionals in which the only purely local term is the ideal 
gas contribution to the free energy, and conversely, this 
theorem is not satisfied by functionals built using a local 
density approximation. 



II. THEORY AND METHODS 

We construct a classical density functional for water, 
which reduces in the homogeneous limit to the opti- 
mal SAFT model for water found by Clark et al. The 
Hclmholtz free energy is constructed using the first four 
terms from Equation [U F id , F hs , F disp and F assoc . In 
the following sections, we will introduce the terms of this 
functional. 
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A. Ideal gas functional 

The first term is the ideal gas free energy functional, 
which is purely local: 



F id [n] = k B T j n(v) fin 



n(r) 



1 dr, (5) 



where n(r) is the density of water molecules and uq is 
the quantum concentration 



mksT 
2irh 2 



3/2 



(6) 



The ideal gas free energy functional on its own satisfies 
the contact value theorem and its limiting case of small 
solutes (Equations [3] and @| . These properties are re- 
tained by our total functional, since all the remaining 
terms are purely nonlocal. 



B. Hard-sphere repulsion 

We treat the hard-sphere repulsive interactions us- 
ing the White Bear version of the Fundamental-Measure 
Theory (FMT) functional for the hard-sphere nuic~&. 
FMT functionals are expressed as the integral of the fun- 
damental measures of a fluid, which provide local mea- 
sures of quantities such as the filling fraction, density of 
spheres touching a given point and mean curvature. The 
hard-sphere excess free energy is written as: 

F hs [n] = k B T J ($i(r) + $ 2 (r) + $ 3 (r))dr , (7) 

with integrands 



$i = -n ln(l - n 3 ) 

nin 2 - nyi • n V2 



1 - 7l 3 



(8) 
(9) 



*3 = («2 - 3n 2 n y2 • n V2 ) 



n 3 + (l-n 3 ) 2 ln(l-?i 3 ) 



36irn 2 (1 — n 3 ) 2 



(10) 



where the fundamental measure densities are given by: 
n 3 (r) = /"n(r')e(|r - r'| - R)dr' (11) 
n 2 (r) = J n(r')8(\r - r'\ - R)dr' (12) 



n V2 = Vn 3 

n 2 



m = 

nyi = 

n = 



4irR 
ny 2 

AttR 
n 2 

4nR 2 ' 



(13) 
(14) 

(15) 

(16) 



The density n 3 is the filling fraction and n 2 describes 
the number of spheres touching a given point. For our 
functional for water, we use the hard-sphere radius of 
3.03420 A, which was found to be optimal by Clark et 

„/34 



C. Dispersion free energy 

The dispersion free energy includes the van der Waals 
attraction and any orientation-independent interactions. 
We use a dispersion term based on the SAFT-VR 
approach^, which has two free parameters (taken from 
Clark et alM.): an interaction energy ed and a length scale 
XdR- 

The SAFT-VR dispersion free energy has the forrn^ 



F, 



disp [ 



(oi(r)+/3a 2 (r))n(r)dr, (17) 



where a\ and a 2 are the first two terms in a high- 
temperature perturbation expansion and /3 = 1/ksT. 
The first term, ai, is the mean-field dispersion interac- 
tion. The second term, a 2 , describes the effect of fluctu- 
ations resulting from compression of the fluid due to the 
dispersion interaction itself, and is approximated using 
the local compressibility approximation (LCA), which 
assumes the energy fluctuation is simply related to the 
compressibility of a hard-sphere reference fluid^. 

The form of a\ and a 2 for SAFT-VR is given in ref- 
erence^, expressed in terms of the filling fraction. In 
order to apply this form to an inhomogeneous density 
distribution, we construct an effective local filling frac- 
tion for dispersion rjd, given by a Gaussian convolution 
of the density: 



Vd(r) 



=-3-3 / n(r')cxp 



„'|2 



2{2\ dSd Rf 



dr'. 



(18) 



This effective filling fraction is used throughout the dis- 
persion functional, and represents a filling fraction aver- 
aged over the effective range of the dispersive interaction. 
Here we have introduced an additional empirical param- 
eter Sd which modifies the length scale over which the 
dispersion interaction is correlated. 



D. Association free energy 

The final attractive energy term is the association 
term, which accounts for hydrogen bonding. Hydrogen 
bonds are modeled as four attractive patches ("associa- 
tion sites" ) on the surface of the hard sphere. These four 
sites represent two protons and two electron lone pairs. 
There is an attractive energy e a when two molecules are 
oriented such that the proton of one overlaps with the 
lone pair of the other. The volume over which this inter- 
action occurs is K a , giving the association term in the free 



4 



energy two empirical parameters that are fit to the ex- 
perimental equation of state of water (again, taken from 
Clark et aM). 

The association functional we use is a modified ver- 
sion of Yu and Wu^, which includes the effects of den- 
sity inhomogeneities in the contact value of the correla- 
tion function g„ s , but is based on the SAFT-HS model, 
rather than the SAFT-VR model^, which is used in the 
optimal SAFT parametrization for water of Clark et aJ^s. 
Adapting Yu and Wu's association free energy to SAFT- 
VR simply involves the addition of a correction term in 
the correlation function (sec Equation |23|). 

The association functional we use is constructed by us- 
ing the density no(r), which is the density of hard spheres 
touching a given point, in the standard SAFT-VR as- 
sociation energy^. The association free energy for our 
four-site model has the form 



Ak B T J n (r)C(r) (lnX{r) 



X(r) 



(19) 



where the factor of 4 comes from the four association 
sites per molecule, the functional X is the fraction of 
association sites not hydrogen-bonded, and £(r) is a di- 
mensionless measure of the density inhomogencity. 



C(r) = f 



n 2 v ■ n 2 y 



a- 



(20) 



The fraction X is determined by the quadratic equation 



A(r) 



y/l + 8n (r)C(r)A(r)-l 
4n (r)C(r)A(r) 



(21) 



where the functional A is a measure of hydrogen-bonding 
probability, given by 



A(r) = n a g s a w {r) (e 



-/3e„ 



1 



S f(r)= S -(r) + ^ 



HSf 



1 / dcii 



X d da-i 



4 \dr) d {r) 3r] d d\ d 



(22) 
(23) 



where g^ w is the correlation function evaluated at con- 
tact for a hard-sphere fluid with a square- well dispersion 
potential, and a\ and a 2 are the two terms in the disper- 
sion free energy. The correlation function g^ w is written 
as a perturbative correction to the hard-sphere fluid cor- 
relation function , for which we use the functional of 
Yu and Wui2 : 



9? S = 



1 



R C,n 2 



R 



(n 2 2 



1 - n 3 2 (f - n 3 ) 2 18 (1 - n 3 )= 



(24) 



E. Determining the empirical parameters 

The majority of the empirical parameters used in our 
functional arc taken from the paper of Clark et al on 
developing an optimal SAFT model for water—. This 
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FIG. 2. Comparison of Surface tension versus temperature for 
theoretical and experimental data. The experimental data is 
taken from NIST.— The length-scaling parameter s d is fit so 
that the theoretical surface tension will match the experimen- 
tal surface tension near room temperature. 



SAFT model contains five empirical parameters: the 
hard-sphere radius, an energy and length scale for the 
dispersion interaction, and an energy and length scale 
for the association interaction. In addition to the five 
empirical parameters of Clark et al, we add a single ad- 
ditional dimensionless parameter s d — with a fitted value 
of 0.353 — which determines the length scale over which 
the density is averaged when computing the dispersion 
free energy and its derivative. We determine this fi- 
nal parameter by fitting the to the experimental sur- 
face tension with the result shown in Figured Because 
the SAFT model of Clark et al overestimates the criti- 
cal temperature — which is a common feature of SAFT- 
based functionals that do not explicitly treat the critical 
point — we cannot reasonably describe the surface tension 
at all temperatures, and choose to fit the surface tension 
at and around room temperature. 



III. RESULTS AND DISCUSSION 

A. One hydrophobic rod 

We begin by studying a single hydrophobic rod im- 
mersed in water. In Figure |3] we show the excess chemi- 
cal potential at room temperature, scaled by the solvent- 
accessible surface area of the hard rod, plotted as a func- 
tion of hard-rod radius. We define the hard-rods radius 
as the radius from which water is excluded. For rods with 
radius larger than 0.5 nm or so, this reaches a maximum 
value of 75 mN/m, which is slightly higher than the bulk 
surface. In the limit of very large rods, this value will 
decrease and approach the bulk surface value. As seen 
in the inset of Fig. [3J for rods with very small radius 
(less than about 0.5 A) the excess chemical potential is 
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FIG. 3. Excess chemical potential per area versus radius for 
a single hydrophobic rod immersed in water. This should have 
an asymptote equal to the surface tension at room tempera- 
ture, and it agrees well with the surface tension in Figure [2] 
The inset for rods with a very small radius shows the linear 
relationship expected based on Equation [4] 
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FIG. 4. Density profiles for single rods of different radii. The 
dotted line represents the saturated liquid density and the 
points represent the expected contact density derived from 
the contact value theorem and calculated free energy data. 



proportional to volume, as required by the contact-value 
theorem (see Equation 2]). 

We show in Figure |4] density profiles for different radii 
rods, as well as the prediction for the contact value of 
the density as a function of rod radius, as computed from 
the free energies plotted in Figure [H The agreement be- 
tween these curves confirms that our functional satisfies 
the contact-value theorem and that our minimization is 
well converged. As expected, as the radius of the rods be- 
comes zero the contact density approaches the bulk den- 
sity, and as the radius becomes large, the contact density 
will approach the vapor density. 



Before transition 




y (nm) 

FIG. 5. Density profiles illustrating the transition from 
vapor to liquid water between the rods. The radius is 0.6 nm, 
the top figure is at a separation of 0.6 nm and the bottom 
is 0.7 nm. Figure [5] shows the energy for these and other 
separations. 

B. Hydrophobic interaction of two rods 

We now look at the more interesting problem of two 
parallel hard rods in water, separated by a distance d, 
as shown in Figure [SJ At small separations there is only 
vapor between the rods, but as the rods are pulled apart, 
the vapor region expands until a critical separation is 
reached at which point liquid water fills the region be- 
tween the rods. Figure [5] shows density profiles before 
and after this transition for rods of radius 0.6 nm. This 
critical separation for the transition to liquid depends on 
the radii of the rods, and is about 0.65 nm for the rods 
shown in Figure [5] The critical separation will be differ- 
ent for a system where there is attraction between the 
rods and water. At small separations, the shape of the 
water around the two rods makes them appear as one 
solid "stadium" -shaped object (a rectangle with semi- 
circles on both ends). 

To understand this critical separation, we consider the 
free energy in the macroscopic limit, which is given by 

F = jA + pV. (25) 

The first term describes the surface energy and the sec- 
ond term is the work needed to create a cavity of volume 
V. Since the pressure term scales with volume, it can be 
neglected relative to the surface term provided the length 
scale is small compared with j/p, which is around 20 /im, 
and is much larger than any of the systems we study. For 
micron-scale rods, the water on the sides of the 'stadium' 
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FIG. 6. Free energy of interaction (also known as the poten- 
tial of mean force) versus separation for two hydrophobic rods 
ranging in radius from 0.2 nm to 1.2 nm. All were arbitrar- 
ily offset to zero at large separations for ease of comparison. 
The transition corresponds to the phase change from vapor 
to liquid between the rods as pictured in the density profiles 
in Figure [5] 

configuration will bow inward between the rods and the 
density will reduce to vapor near the center point where 
the rods are closest to each other. 

Starting from the surface energy term, we can calculate 
the free energy per length, which is equal to the circum- 
ference multiplied by the surface tension. The force per 
length is the derivative of this with respect to the sepa- 
ration. The circumference of the stadium-shape is 

C s = 2nr + 4r + 2d (26) 

and so the force per length is equal to twice the surface 
tension. 

We plot in Figure |6] the computed free energy of in- 
teraction per unit length from our classical density func- 
tional (solid lines) , as a function of the separation d, along 
with the free energy predicted by our simple macroscopic 
model (dashed lines). The models agree very well on the 
force between the two rods at close separations, and have 
reasonable agreement as to the critical separation for rods 
greater than 0.5 nm in radius. 

Walther et a&, studied the interactions between two 
carbon nanotubes, which are geometrically similar to our 
hydrophobic rods, using molecular dynamics with the 
SPC model for water, and a Lennard- Jones potential for 
the interaction of carbon with water, for nanotubes of 
diameter 1.25 nm and separations ranging from about 
0.3 nm to 1.5 nm. The SPC model underestimates the 
surface tension of water by about 24%^, so we cannot 
expect this work to provide quantitative agreement with 
real water. Walther et al observe a drying transition 
between the two nanotubes, which occurs at a smaller 
radius than our results suggest. However, when Walther 
et al disable the attractive interaction between nanotubc 



FIG. 7. Free energy of interaction versus separation for four 
hydrophobic rods ranging in radius from 0.2 nm to 1.2 nm. 
All were arbitrarily offset to zero at large separations. The 
transition corresponds to the phase change from vapor to liq- 
uid between the rods as pictured in the density profiles in 
Figure [8] 

and water, the drying effect occurs at much longer range, 
in agreement with our results. 

C. Hydrophobic interactions of four rods 

We go on to study four parallel hard rods, as examined 
by Lum, Chandler and Weeks in their classic paper on 
hydrophobicity at small and large length scales^. As in 
the case of two rods — and as predicted by Lum et al — 
we observe a drying transition, as seen the density plot 
shown in Figure [8] In Figure El we plot the free energy 
of interaction together with the macroscopic approxima- 
tion, and find good agreement for rods larger than 0.5 nm 
in radius. This free energy plot is qualitatively similar to 
that predicted by the LCW theory^, with the difference 
that we find no significant barrier to the association of 
four rods. 



D. Hydration energy of hard-sphere solutes 

A common model of hydrophobic solutes is the hard- 
sphere solute, which is the simplest possible solute, and 
serves as a test case for understanding of hydrophobic 
solutes in water—. As in the single rod, we begin by 
examining the ratio of the excess chemical potential of 
the cavity system to the solvent-accessible surface area 
(Figure [5]) . This effective surface tension surpasses the 
bulk surface tension at a radius of almost 1 nm, and at 
large radius will drop to the bulk value. As with the sin- 
gle rod, we see the analytically correct behavior in the 
limit of small solutes. For comparison, we plot the free 
energy calculated using a molecular dynamics simulation 
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FIG. 8. Density profiles illustrating the transition from 
vapor to liquid water between four rods. The radius is 0.6 nm, 
the top figure is at a separation of 1.53 nm and the bottom 
is 1.56 nm. Figure [7] shows the energy for these and other 
separations. 



of SPC/E water—. The agreement is quite good, apart 
from the issue that the SPC/E model for water signifi- 
cantly underestimates the surface bulk tension of water 
at room temperature^. 

Figure [10] shows the density profile for several hard 
sphere radii, plotted together with the results of the 
same SPC/E molecular dynamics simulation shown in 
Figure [Op 4 . The agreement with simulation is quite rea- 
sonable. The largest disagreement involves the density 
at contact, which according to the contact value theorem 
cannot agree, since the free energies do not agree. 
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FIG. 9. Excess chemical potential per area versus radius 
for a single hydrophobic sphere immersed in water. This 
should have an asymptote equal to the surface tension at 
room temperature, and it agrees well with the surface ten- 
sion in Figure [2] Results from a simulation of SPC/E watep^ 
are shown as circles. The horizontal lines show the experi- 
mental and SPC/E bulk surface tension for water at standard 
atmospheric temperature and pressure. 
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FIG. 10. Density profiles around hard-sphere solutes of dif- 
ferent radii. Predictions from our classical density-functional 
theory are in solid red, while the dotted line shows the result 
of a molecular dynamics simulation of SPC/E watei^. 



IV. CONCLUSION 

We have developed a classical density functional 
for water that combines SAFT with the fundamental- 
measure theory for hard spheres, using one additional 
empirical parameter beyond those in the SAFT equation 
of state, which is used to match the experimental surface 
tension. This functional does not make a local density 
approximation, and therefore correctly models water at 



both small and large length scales. In addition, like all 
FMT functionals, this functional is expressed entirely in 
terms of convolutions of the density, which makes it effi- 
cient to compute and minimize. 

We apply this functional to the case of hard hydropho- 
bic rods and spheres in water. For systems of two or 
four hydrophobic rods surrounded by water, we see a 
transition from a vapor-filled state a liquid-filled state. 
A simple model treatment for the critical separation for 
this transition works well for rods with diameters larger 
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than 1 nm. In the case of spherical solutes, we find good agreement with SPC/E simulations. 
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